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Abstract 



A test particle code is used to study the response of ensembles of particles to a 
two-dimentional, time-dependent model of the geomagnetic tail, and test the proposition 
(Coroniti, 1985a, b; Buchner and Zelenyi, 1986; Chen and Palmadesso, 1986; Martin, 1986) 
that the stochasticity of the particle orbits in these fields is an important part of the 
physical mechanism for magnetospheric substorms. The realistic results obtained for the 
fluid moments of the particle distribution with this simple model, and their insensitivity 
to initial conditions, is consistent with this hypothesis. 
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I Introduction 


Since Dungey (1961) posed the problem of magnetic field line reconnection in the 
Earth’s magnetotail a search has been on for a mechanism that would allow the release of 
the stored magnetic energy, leading to a magnetospheric substorm. The tearing mode is 
a popular candidate, but since the magnetotail is highly collisionless, we cannot count on 
collisions to provide the dissipation necessary for the instability, and the reconnection time 
given by resistive MHD theory (Furth et. a il. 1963) is several orders of magnitude too long. 
Therefore, either a collisionless resistivity is required or a mechanism of an ideal MHD type 
for reconnection, such as the coalescence instability. We focus on the former in the present 
paper. In the search for that resistivity two broad approaches can be identified: a search 
for a dynamical collisionless resistivity and a search for a turbulent collisionless resistivity. 
In the former, the finite coherence time (the time during which the current carriers are 
coherently accelerated by the electric field, also called the decorrelation time) arises from 
single particle dynamics in the magnetotail fields, while in the latter it comes from the 
scattering of the particles by the modes of the current carrying plasma. In the present 
work we consentrate on dynamical collisionless resistivity. Several mechanisms have been 
proposed for turbulent collisionless resistivity (eg. Hagege et. al., 1973; Huba et. al., 1978; 
Drake and Lee, 1977; Esarey and Molvig, 1987), but recent results form the ISEE-1 and 
ISEE-2 satellites on the CDAW-6 substorm (Anderson, 1984; Dusenbery, 1989) indicate 
that wave intensities in the current sheet are insufficient to account for the observed recon- 
nection rate and, in addition, most wave intensities decrease dramatically at the center of 
the current sheet, making it even harder for turbulence to produce the required resistivity. 

Measurements of the particle distribution function by the ISEE-1, ISEE-2 and ISEE-3 
satellites indicate the presence of unmagnetized ion orbits in the plasma sheet and the 
earthward flowing ion streams in the plasma sheet boundary layer (PSBL). On the basis 
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of these data, Lyons and Nishida (1988) recently suggested that the onset of the neutral 
line earthward of 80i2 e involves the details of this diffusion zone in the plasma sheet where 
the ion motion is non-adiabatic. Their suggestion is that the neutral line must form in this 
ion diffusion-like zone and that this plasma is also the source of the observed earthward 
streaming ions. In addition, Chen and Palmadesso (1986), Martin (1986) and Buchner and 
Zelenyi (1986), recently demonstrated using a variety of techniques (Poincare surface of 
section plots, Lyapunov exponents, separatrix crossing) the stochastic nature of the particle 
orbits in the parabolic field lines of the geomagnetic tail. As is shown in Part- II below, 
the stochasticity arises from the resonance overlap of the nonlinear current sheet cyclotron 
oscillations in the reversed tail-earthward field with the cyclotron motion in the weak 
northward dipole field component B z . Coroniti (1985a) and Buchner and Zelenyi (1987) 
also proposed substorm models that rely entirely on single particle dynamics. 

In view of these new developments, we have undertaken a study of the dynamics and 
transport of particles in the parabolic field lines of the geomagnetic tail in the presence of 
a neutral line, at a distance of 25-35 R e , produced by a reconnection mode perturbation. 
We advance ensembles of particles with a test particle code, so that no collective effects 
are present other, than the ones implicit in the assumption that the overall MHD picture 
is correct (eg. Steinolfson and Van Hoven, 1984) and the perturbation is given by a re- 
connection mode. All the phenomena we observe are the result solely of the dynamics of 
ensembles of charged particles in the model fields. The electric and magnetic fields in the 
geomagnetic tail are modeled by 


B x — Bo tanh^-^ 

B z = B 0 b 0 + B 0 tjjoksin(kx)e li (1) 

E y = Bo^o- cos(fcx)e 7 * 
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with x,y, z the usual magnetotail coordinates. These fields can be derived from the mag- 
netic flux function 


ip — Bo [a ln(cosh(z/a)) + V’o cos(fcx)e 7t — 60 ^] ( 2 ) 

where the vector potential is A = -rp(x, 2 , t)e y and the fields are given by B = y x Vxp and 
E y = (l/c)d t tp. In the above expressions B 0 is the magnitude of the asymptotic (far from 
the neutral sheet) horizontal field, a the scale length of the horizontal field, &o the normal 
(northward) magnetic field component, V’o I s the amplitude of the tearing perturbation, 
7 its growth rate and k its wavenumber. The electromagnetic perturbation in Eq. (2) 
is the standard tearing mode perturbation well studied in resisitive MHD simulations 
(eg. Steinolfson and Van Hoven, 1984). The neutral line occurs when B z = -dtp/dx = 0, 
i.e. at time t c 



for the model fields in Eq. (1). For t < t c the field lines are long, open loops. 

The particles are advanced in time using the fields given in Eq. (1), and particle 
density, bulk velocities, temperature and current are measured as a function of x and 
2 . We find that the initial distribution changes little until the neutral line is formed, 
and then the number density in the vicinity of the X-point drops, a narrow channel of 
cross-tail current, in phase with the electric field, is formed around the X-point, and the 
temperature increases in the same narrow channel where the cross-tail current flows, with 
final temperatures in the range of 1-12 keV. Finally the hot particles are expelled from the 
vicinity of the X-point into the PSBL, and earthward streaming ions with bulk velocities 
in the range of 400-1000 km/s are observed there. These results are in good qualitative 
and quantitative agreement with observed events in the magnetotail. Number density is 
observed to decrease with a simultaneous increase in temperature after substorm onset, 
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with temperatures in the range of 1-10 keV (eg. Paschmann et. al., 1985; Huang, 1989), 
and earthward streaming ions with velocities of several hundred km/s are observed in the 
PSBL (eg. DeCoster and Prank, 1979; Eastman et. al., 1984). The agreement of our results, 
which were obtained with a test particle code, with the observations lends support to the 
proposition (Coroniti, 1985a; Buchner and Zelenyi, 1986; Martin, 1986) that single particle 
dynamics play an important role in determining the behavior of the geomagnetic tail. 

II Previous Work 

Harris (1962) solved the Vlasov equation for the case of a reversing one-dimensional 
magnetic field, a configuration often used as a zeroth order approximation to the mag- 
netotail fields. More realistic approximations include a cross-tail electric field and/or a 
small normal magnetic field component, and many authors have studied the dynamics of 
single particles in these model fields (eg. Parker, 1957; Speiser, 1965, 1967, 1968, 1970; 
Schindler, 1965; Sonnerup, 1971; Eastwood, 1972, 1974; Lyons and Speiser, 1982; Kim 
and Cary, 1983; Speiser and Lyons, 1984; Martin, 1986; Buchner and Zelenyi, 1986, 1987). 
The first attempt to calculate a collisionless resistivity, and hence the tearing mode growth 
rate, was made by Coppi et. al. (1966) and Laval et. al. (1965) for a one-dimensional Harris 
sheet. They found the system unstable to the electron tearing mode, but the growth time 
was unrealistically long for magnetospheric parameters, ~ 4.5 days. The more realistic case 
of a two-dimensional magnetic field (a small B z added) was first studied by Galeev and Ze- 
lenyi (1976) ( see also Galeev, 1979) who found that, as was proposed by Schindler (1974), 
the electron tearing mode is stable in this field geometry. This is due to the fact that the 
normal magnetic field component magnetizes the electrons in the quasineutral layer (where 
the horizontal field vanishes) and forces an adiabatic response to the tearing perturbation. 
The electron Landau interaction, which provides the dissipation necessary for the growth 
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of the tearing mode in the one- dimensional case, is inhibited and the mode is therefore sta- 
ble. From the energy balance point of view, the energy required to produce the adiabatic 
electron response to the tearing mode exceeds the free energy available in the current sheet 
(Coroniti, 1980; Lembege and Pellat, 1982). Furthermore, the calculation of Lembege and 
Pellat (1982), which included the effects of the particle drift in the normal field as well as 
a finite Larmor radius, showed that the adiabaticity of the electrons will even stabilize the 
ion tearing mode, which was proposed by Schindler (1974) as an alternative to the electron 
tearing mode, and the growth rate of which was calculated by Galeev and Zelenyi (1976). 
It was therefore becoming important to search for a mechanism that could destabilize the 
electron tearing mode. 

As was mentioned earlier this search can be divided into two broad categories; dy- 
namical or inertial resistivity and turbulent or collective resistivity. Coroniti (1980) first 
considered the effects of turbulence-induced pitch-angle scattering of the electrons on the 
growth rate of the electron tearing mode, and found that pitch- angle diffusion makes the 
mode unstable. Many authors (eg. Hagege et. al., 1973; Huba et. al., 1978; Drake and 
Lee, 1977; Esaxey and Molvig, 1987) have considered the effects of different types of turbu- 
lence (lower hybrid drift wave turbulence, ion acoustic turbulence, broadband electrostatic 
noise) and they all cause the electron tearing mode to become unstable. Turbulence is 
indeed present in the geomagnetic tail, but it is either in the wrong place (broadband elec- 
trostatic noise, for example, is mostly confined to the PSBL) or of insufficient magnitude. 
In addition, recent studies (Anderson, 1984; Dusenbery, 1989) show that wave intesity 
decreases even more at the center of the current sheet. This situation, together with the 
fact that all kinds of random driving considered (including spatial rather than velocity- 
space diffusion) destabilized the electron tearing mode, led to the suggestion (Buchner and 
Zelenyi, 1986; Martin, 1986; Doxas, 1988) that the intrinsic stochasticity of the particle 
orbits in the magnetotail fields may be sufficient to destabilize the mode. 
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Kim and Cary (1983) first showed that particle trajectories in the simpler but related 
geometry of elliptic magnetic flux surfaces can be stochastic. That was also demonstrated 
for more realistic magnetotail field models by the work of Chen and Palmadesso (1986) 
who found stochastic regions in phase space. Martin (1986) measured Lyapunov exponents 
and used them as the effective coherence time r to calculate the collisionless resistivity. 
Buchner and Zelenyi (1986, 1987) identified 

K = ^»/Z=^ (4) 

B 0 V Po \A 

as the parameter that controls the stochasticity of the orbits. In the above expression B z q 
is the magnitude of the constant normal magnetic field component, Bo the asymptotic (far 
from the neutral sheet) value of the horizontal magnetic field, a the scale length of the 
reversing horizontal field (which is equal to the width of the Harris current sheet) and po 
the particle gyroradius in the asymptotic field J?o- The parameter e = po/ a is the small 
parameter in the adiabatic approximation. In the current sheet the parameter k arises 
from the ratio of the cyclotron frequency u> cz = eB 2 o/mc around the north-south field, 
to the frequency of the strongly nonlinear vertical oscillations in the reversed earth-tail 
magnetic field u> C x( z ) ~ cB x {z)/mc (eg. Lyons, 1984). The value of k gives the number 
of gyrations around the normal field that a particle will execute in one period of the 
vertical oscillations. For k 1 the magnetic moment p is a good adiabatic invariant, 
and the adiabatic approximation is valid. For k <C 1 a new adiabatic invariant can be 
found (Schindler, 1965; Speiser, 1970; Sonnerup, 1971). For k « 1 however, no adiabatic 
invariant exists, and the motion is highly stochastic. 

Ill Numerical Results 

The equations of motion used to advance the particles are given by the Lorentz force 
with the fields in Eq. (1). We introduce dimensionless variables by scaling length to the 
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perturbation wavelength A = 2ir/k, time to the ion gyroperiod in the asymptotic field To,- = 
27r/u;oi = 27r(mic/ei?o) = 3.277 sec, mass to the ion mass and charge to the unit charge e 
(electrons have charge — e). The unit velocity is therefore v 0 = &oi/k = 9.71 x 10 9 cm/s , 
and the unit current jo = ep n v o = 1.55 pA/rn 2 , where p n — 0.1 cm~ 3 is the number 
density in the magnetotail, and a typical value for the asymptotic field is Bo = 20 nT 
(eg. Lyons and Williams, 1984). From now on till values given for the parameters and all 
variables in equations will be dimensionless unless otherwise stated. Whenever needed to 
avoid confusion, the dimensionless quantities will be denoted by a hat (eg. x = x/A). The 
system is described by four parameters. The amplitude of the magnetic flux associated 
with the tearing mode tpo, the constant normal magnetic field component bo, the aspect 
ratio of the tearing mode ka(= 27ra) and its growth rate 7. 

We fully resolve the smallest time and space scales in the problem. For the parameters 
used, the smallest length scale is the particle gyroradius in the asymptotic field, po = v/u>o 
with v the total particle velocity. The length scale over which the magnetic field reverses 
direction, a, is always at least two orders of magnitude larger than p even in very high 
energy test runs, and the wavelength of the tearing perturbation is the longest scale in 
the system. Approximate values for these parameters in the magnetotail are (eg. Lyons 
and Williams, 1984) A = 8-80 R e and a fa 1 R e where R e = 6380 km is the radius of the 
Earth. In contrast the electron gyroradius is of the order of a few kilometers, and the ion 
gyroradius a few hundred kilometers. The smallest timescale is the particle gyroperiod 
in the asymptotic field To = 27r/a>o = 1. The only other timescale in the problem is the 
tearing mode growth time, and the shortest growth time we used was 20 ion gyroperiods. 

Sample populations of 16000-64000 particles are advanced in time using a fourth order 
Runge-Kutta method (eg. W. Gear, 1971; W.H. Press et. al., 1986). Electrons and ions 
are treated separately so that a realistic ion to electron mass ratio of 1836 can be used. 
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The electron gyrofrequency is therefore u>oe = —1836 x (27r). The time step h is chosen so 
that the particle gyroperiod is adequately resolved and the final (cumulative) error in the 
particle position \8x\ is much smaller than the particle gyroradius in the asymptotic field. 
To get an estimate for the error in the particle position, a small percentage of particles 
1%) are also advanced with a sixth order Runge-Kutta and half the timestep. The 
particle position calculated with the higher order method and shorter timestep is 

then compared to the position of the particle Xh calculated with the regular method and 
timestep, and the error \8x\ is given by 

|(5x| = \x h /2 - x h \ (5) 

The lowest value of To/h used in any run is 110, and the lowest value of p/\8x\ is of the 
order of 10 2 . The following values of the parameters were used in the runs: the normal 
magnetic field in the range 0.003 < b Q < 0.1, the amplitude of the tearing mode in the 
range 4 x 10“ 4 < V>o < 0.01, the tearing mode growth rate in the range 0.005 < 7 < 0.05 
and the tearing mode aspect ratio was A/a = 50. In order to test heating ireversibility, as 
well as adress some concerns explained below, we also made a number of runs with a time 
dependent perturbation amplitude (i/>o — ► ^oC Qi and t/>o ^o s hi(u>f)) with —0.001 < ol < 
0.01 and 0 < uit < 71*. 

The equations of motion are periodic in the x-direction so the particles are kept in 
a cell with an x-dimension of one, but are allowed to move freely in z. The particles are 
loaded at t = 0 using two different initial distributions; a flat distribution in both velocity 
and configuration space, 

— X < X < X* , — Z* < Z < Z 9 , — V 1 < Vj < v* (^ a ) 


or a gaussian distribution in z and v, 

1 


p(z,v) 


r-z 2 /2(T Z -v 2 /2<r v 

*.* 5 ( 2 * ) 2 


— x' < x < x 1 


( 66 ) 
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We always start with a flat distribution in x since the equation of motion are periodic in x. 
The range of initial conditions used in the different computer experiments are x 1 = A/2, z 
and <r z in the range 0.1a ^ z 1 , o z 2a, and the velocity space limit v 1 and deviation o v are 
chosen so that the average kinetic energy of the particles is in the range 2 x 10 -5 keV <T< 
2 keV. The system is typically advanced for 200-1000 gyroperiods, which is equivalent to 
10-50 minutes for ions and 0.5-2 seconds for electrons. Time enters explicitly the equations 
of motion exclusively through the combination V’oe 7 *, (we always consider a constant t/>o, 
the time dependence is introduced explicitly later) so that a sweep in the value of xp 0 is 
equivalent to running longer times, with only the initial conditions distinguishing the two 
cases. This observation is supported by the simulations and was used in several occasions 
to shorten lengthy runs. Electron runs were spaced at various times in the 10-50 minute 
time history of the ions using different values of xj) o, so that ‘snapshots of the electron 
behaviour during the ion run were obtained. 

We compute and study the average density, velocity (x, y and z-component), current 
(x, y and z-component), temperature and (j ■ E) as a function of x and z, at any given 
time t. 

The average velocity is defined by 

(Vj)(s,t) = = v 0 (vj)(s,t) (7) 

where n is the number of particles whose s-coordinate is in the range s — > s + 6s irrespective 
of the value of the other coordinate, and i is the j-component of the dimensionless velocity 
of the i th particle in units of the scale velocity vo = u>o i/k. Here s represents either x or z. 

Using the average velocity, we can define the current by 

jj(s,t) = ep p (s)(vj)(s,t) (8) 
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where { Vj ) is the average velocity of all particles whose ^-coordinate has a value between 
s and s + 8$, and p p (s) is the particle number density in the cell s -+ s -{- 8s. The number 
density is in turn given by 

ft-W = (9) 

where n has the same meaning as in Eq. (7), N is the total number of particles in the 
simulation, p n = 0.1cm -3 the average number density in the magnetotail and Ss the 
dimensionless cell size. The factor n/NSs is a density factor that describes how the density 
at position s differs from the average magnetotail number density p n . If n = NSs, then 
p p = p n . With that definition of the density we can now rewrite Eq. (8) as 


r^n ~ i 

ij(M) = ep p (vj)(s,t) = 1 

= epnV ° (jkl ^ = 


( 10 ) 


where jo = ep n v o is the scale current. 

The in-phase part of the reconnection current (j • E) is computed by 


O' ■ £>(..!) = ( e Pn V ^ 5 i COS ( 2,ri ’)^ 


= (epnVo)(j • E) 


( 11 ) 


and the temperature is given in keV using 


im(u 2 ) = ^mvl(v 2 ) 


< 


V 2 ) = {[v x - (6*» 2 ^ + ({?y - {Vy)f^ + ^(”z - ^z)) 2 ) 


(12) 


We have fifteen ion runs and four of the much more expensive electron runs. Although 
nineteen runs are insufficient for a systematic exploration of the parameter space (there 
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are four parameters that describe the system, ifio, ^o, ka rind 7), their qualitative results 
are very similar even for widely different values of the parameters and initial conditions, 
and we feel that we can draw broad conclusions about the behavior of the system even 
from a limited sample. We will go here in detail through two runs (one for ions and one 
for electrons) that exhibit the main behavior patterns of the system. 

For the ion run, the dimensionless variables are t/’o = 0.001, bo = 0.05, 7 = 0.01, 
and a = 0.02. Substituting these values into Eq. (3), we find that the neutral line will 
appear at time t = 207.4 (the dimensionless time is in units of ion gyroperiods, To, = 
2n/u>oi = 3.28 sec). In Fig. (1) we show plots of the magnetic field lines, with a sample 
200 particles superimposed, for three different times during the run. The dashed vertical 
lines are located at z = ±a. As we see in Fig. (l)-a, the particles are loaded at time t = 0 
with a flat distribution between z = —a and z = a, (z* = a in Eq. (6)) and over a whole 
wavelength in x (x f = A/2). The initial velocity distribution is also flat in all velocity 
components (— v'j < vj < Vj) with v) = 1.4 x 10~ 4 , equivalent to an initial temperature of 
only 0.8 eV. As we shall see, even this run, with a very low initial temperature, will give 
results similar to runs with initial temperature close to the observed temperature. 

The initial distribution changes little until shortly before the neutral line forms. In 
Fig. (l)-b we show a plot of the magnetic field lines with 200 sample particles superimposed 
at time t = 183.3, less than 25 ion gyroperiods before the neutral line appears. We see 
that the area around x = —A/4, where the neutral line is first going to form, is thinning 
out, while all particles have been expelled from the centers of the magnetic islands at the 
top and bottom of the picture. In Fig. (l)-c we see the same area at time t = 233.3, only 
26 gyroperiods after the appearance of the neutral line. We see that the area around the 
neutral line is almost depleted. Still, most particles are within the z = ±a boundaries 
they were placed in at t = 0. Fig. ( 1 )-d shows the fields and particles at t = 280.0. All 
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particles have been expelled from the neutral line area, and they are following the magnetic 
field lines, as the magnetic island expands, into the PSBL. Our run is terminated shortly 
afterwards, since we can no longer make any meaningful measurements of the current or 
temperature in the neutral line area. It is interesting to note how remarkably little happens 
in the 180 or so gyroperiods (~ 10 min) at the beginning of the run, while most changes 
take place within 50-70 gyroperiods 3 min) around the appearance of the neutral line. 

The in-phase current (j • E) and the particle temperature is plotted in Figure (2) 
as a function of z for three different times. The initial distribution is again essentially 
unchanged for the beginning of the run. As we see in Figs. (2)-a,b, at time t = 166.7 there 
is a well defined current channel in z, and the temperature has increased by a factor of 
3-4 in the same narrow channel in which the current flows. Notice that the temperature 
profile in z has just developed a new narrow channel around z = 0, where the temperature 
is going to increase dramatically between now and the appearance of the neutral line. 
Figs. (2)-c,d show the current and temperature at time t = 216.7, just nine gyroperiods 
after the appearance of the neutral line. The temperature in the narrow channel around 
the reversal layer has increased by almost an order of magnitude, while the rest of the 
temperature profile in z has broadened even more. Finally in Figs. (2)-e,f we show the 
current and temperature profiles at time t = 273.3. Although some current is still flowing, 
it will soon disappear as the last few remaining particles are being expelled from the area. 
The temperature profile in z has already crashed in the center, with all the energetic 
particles expelled from the center and injected into the PSBL where the temperature is 
now a few keV. Figs. (2)-g,h show the in-phase current (j • E) and the particle temperature 
at t = 166.7 plotted against x. The current profile has developed a sharp gradient at the 
point where the neutral line is going to appear (at x = —A/4) and the temperature profile 
follows it closely, with high temperatures confined to earthward of the neutral line. Both 
the magnitude of the in-phase current and the temperature increase after the appearance 
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of the neutral line, but the profiles remain essencially unchanged until the particles are 
expelled from the neutral line region. 


Although we start with an extremely low temperature (0.8 eV), we end up with tem- 
peratures of a few keV, comparable to the actual temperatures in the tail. When we start 
with more realistic temperatures, (eg. 1.0 keV) the overall heating is much lower than the 
four orders of magnitude we saw in this run, and the final temperature is a little higher 


(8-10 keV ) but still close to the observed values. If we start the run at a time closer 
|time the neutral line appears (by using a larger ipo) the system still ends up with 


mately the same temperature. This tendency of the system to end up at essentially 
le state around the time of the appearance of the neutral line is general and not 
l to the temperature, as we already noticed with the density. 

ig. (3) we show the time history of the x (earthward-tailward) component of the 
elocity. At t = 33.3, there is a small (approximately 15 km/ sec) earthward drift 
ow channel around the reversal layer at z = 0 (Fig. (3)-a). At later times, both 
ig. (3)-b) and shortly after the appearance of the neutral line (Fig. (3)-c), the 


ets a little wider and the velocity is higher, but still quite low, of the order of 40- 
The narrow depression around z — 0 is present in almost all our simulations, 
bably due to the fact that the particles follow so-called Speiser orbits. As they 
and B z (Speiser, 1965) their average ar-velocity is zero (f* v± cos Odd = 0). The 
does not go to zero in our plots because we sample over a finite width in z. 
articles are expelled from the neutral line region and into the PSBL (Fig. (3)-d) 
ard streaming velocity rises quickly to approximately 900 km /sec. Notice how 
ard velocity of the particles stays almost unchanged for the 200 To, before the 
appears, and then sudenly, in only 80 To,-, it increases by more than an order 
ude. Again, as in the case of the particle temperature, the qualitative picture is 


insensitive to initial conditions, and again most of the change comes after the neutral line 
has formed. 

Finally, for the ion run, we show in Fig. (4) contour plots of the particle distribution 
function for the five boxed areas of particle concentration in Fig. (l)-c. Figs. (4)-a,-d 
correspond to the four ‘arms’ of particle concetration clockwise from the top left (—1.5 < 
z/a < 0., 0.2 < x/\ < 0.5, etc.). Fig. (4)-e corresponds to the particle concetration 
earthward from the neutral line at —0.5 < z/a < 0.5, —0.2 < x/X < 0.2. The figure 
demonstrates once again the importance of the neutral line, since stochastic velocity-space 
diffusion is seen to occur only in its vicinity. Since the particle distribution function is 
broadenning, the figure also provides evidence that the tremendous rise in temperature we 
observe is real heating, and not the artifact of some unusual distribution function, or a 
reversible flow of energy from the fields to the particles. The question of reversibility is an 
important one, and will be addressed in more detail later. 

The above scenarios for the changes in the density, bulk velocity and temperature are 
in good agreement with observations. Earthward streaming ions with speeds up to 40 km/s 
have been observed in the plasma sheet shortly after substorm onset (eg. Paschmann 
et. al., 1985), while earthward streaming ions with speeds of 100-900 km/s have been 
observed in the PSBL (eg. DeCoster and Frank, 1979; Eastman et. al., 1984). The ion 
temperature is also seen to increase with a simultaneous decrease in density shortly after 
onset (eg. Paschmann et. al., 1985; Huang, 1989) and ion temperatures in the plasma sheet 
are in the range 1-10 keV, although the temperatures observed in the PSBL (« 1 keV ) are 
lower than those seen in the simulations. 

For the electron run the dimensionless parameters are b 0 — 0.003, a = 0.02 and 
7 = 0.01. Following an electron run over several ion gyroperiods is expensive, so we use 
the fact that time and tpo are equivalent to ‘take snapshots’ of the system at different times. 
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This particular run was divided into two parts. In the first we use rpo = 4.0 x 10 — 4 , which 
places the system 17.7 To, before the appearance of the neutral line and in the second we 
use %l > o = 6.0 x 10 -4 , which places the system 22.8 Toi after the appearance of the neutral 
line. We will henceforth disregard the value of ipo, and refer to the two parts of the run 
by the equivalent time label. The run at t = —17.7 for instance, or the early run, is the 
part with V’o = 4.0 x 10~ 4 . Similarly the late run is the part with V’o = 6.0 x 10 -4 . In 
both cases the initial distribution was flat in both configuration and velocity space, with 
v'j = 1.4 x 10 -4 , which is equivalent to an initial temperature of 5 x 10“ 5 keV. 

Both parts of the electron run are similar with the late (post neutral line) run develop- 
ing quicker, and having consistently a much higher current and temperature, as expected. 
The main difference in the two runs, which shows most dramatically the effect of the neu- 
tral fine, is the time history of the in-phase current (j ■ E). For the early electron run, 
the value of (j • E) is plotted against time in Fig. (5)-a, and for the late run in Fig. (5)-b. 
Because of the similarity of the plots in the beginning, it it fair to assume that during that 
time the system is reaching its ‘preferred state’, while after the initial ‘equilibration’ we 
can see the real response of the system. In the early run (j • E) reaches a constant value 
(except for small oscillations), while in the late run it increases almost linearly in time. 
Since the only difference in the two runs is the existence of the neutral line in the late run, 
we see that the appearance of the neutral line causes the value of the height-integrated 
(j ■ E) to increase rapidly with time. Again the picture is insensitive to initial conditions. 
In addition, an otherwise identical pair of runs with a modest amount of velocity-space 
diffusion ( D v = 0.6 km 2 /s 3 ) added in the form of random ‘kicks’ in the particle velocity, 
exhibited virtually identical behavior, as can be seen in Figs. (5)-c,d. A number of runs was 
made with different initial conditions, both flat and gaussian, and velocity-space diffusion 
coefficients as high as 1.6 km 2 /s 3 , and they all showed the same basic qualitative features 
with the corresponding run that had no velocity space diffusion. 
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We therefore see that single particle dynamics can give realystic results for the mag- 
netotail evolution, and that the neutral line plays an important role in this evolution. Two 
important questions remain however. 

1) Is the particle energization truly irreversible, or are we observing some ‘sloshing’ mode, 
and the particles will give up their energy back to the fields once the direction of the 
electric field is reversed, and 

2) How will collective effects change our results? In particular, will the waves become 

depleted as their energy is imparted to the particles and the growth slow down, or will, 
conversly, the currents produced by the particles reinforce the waves? This question be- 
comes particularly important when it is observed that the background current density is 
(c/47t)V 2.5 x 10 -3 \iAfm 2 , much smaller that the currents produced by the ions 

(cf. Fig. (2), the scale current is jo = 1.55 fiA/m 2 ) although still much higher than those 
produced by the electrons (cf. fig. (5)). In addition, the current is now flowing over a much 
narrower channel in z (cf. Fig. 2). When the difference in the width of the current channels 
is taken into account, the total ion current is still comparable to the total MHD current, 
but no longer dominant. 

The first question is partly addressed by the uj_~U|| plots in Fig. (4). The second 
question can of course be fully addressed only in the context of a fully self consistent 
simulation. However, the method we use to show the irreversibility of particle heating, 
can also be employed to provide some useful, albeit imperfect, insight into the question of 
collective effects. 

In order to show that the particle heating is irreversible, we introduce an explicit time 
dependence in the amplitude of the tearing mode. We do this in two different ways, and 
the results we obtain are identical for both methods. In the first method we put 

V>o — * i/>o sin(u><) (13a) 
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in Eq. (1), which, if we also choose 7 = 0, gives 


B z — Bobo + Boipok sin(fcx) sin(u>f) 
u cos (u>t) 


Ey — B 0 iPo- 


cos(kx) 


(136) 


If we now choose u> and ipo appropriately and stop our simulation at ut = 7r/2, the electric 
field has changed direction half way through the run (before or after the appearance of 
the neutral line, depending on the value of ip o), while the magnetic field has maintained 
the same direction throughout the run (the magnetic island increases at first and then 
decreases, but the field topology remains unchanged). In the second method, we put 


ipo V’oe 


af 


(14a) 


in Eq. (1), which gives 


B z = B 0 bo + B 0 iPoksm(kx)e Qt2+ ' llt 
(2 at -f 7) 


(146) 


E y = Boipo' 


cos (kx) 


By appropriately choosing the values of a and 7, we can arange for the electric field 
to change direction (at t r — —7/20) before or after the appearance of the neutral line 
(at t c = (I/7) ln(6 0 /A:^o)) cf Eq. 3) without changing the topology of the magnetic field. 
Several runs were made with these two different methods, using different initial conditions, 
and the results confirm that the heating is indeed irreversible. Although the heating rate 
was substantially reduced, and even stopped, after the electric field changed direction, 
the particles retained the temperatures they had already attained, demonstrating that 
stochasticity can indeed lead to thermalization of the particle energy, rather than just 
streaming. 

The question of depletion or reinforcement of the wave by the particle currents is 
equivalent to the question of what is the proper growth rate of the tearing mode, since 
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there is no fundamental difference between a time dependent amplitude, and a different 
mode growth rate. Any algebraic time dependence will, of cource, be quickly masked by 
the exponential. A simple exponential dependence is equivalent to running different values 
of 7, since we can put 7 = 71 +72 and adopt the view that the wave is growing as e 7l< while 
its amplitude changes as e 75 *. We used a number of different values for 7, in the range 
0.005 < 7 < 0.05. The heating rate was slightly lower for lower 7’s, the fined temperatures 
were also lower (eg. ~ 1 keV for 7 = 0.005 and 3.8 keV for 7 = 0.05 for two otherwise 
identical runs), and the timescales were different as expected (we had to integrate for longer 
times to obtain similar results), but the overall picture remained unchanged. In particular 
the relative importance of the neutral line remained the same for all values of 7. 

An intuitively more obvious, although still rudimentary, way of addressing the question 
is to use the second method mentioned above for the irreversibility of heating, by viewing 
the e"* 2 time dependence as originating from the amplitude of the tearing mode (Eq. 14a). 
A negative value of a would then correspond to a wave being depleted by transfering 
its energy to the particles, while a positive a would correspond to the particle currents 
reinforcing the wave. Several values of a were used in the range —0.001 < a < 0.01, and 
the results were similar to the ones obtained by running different values of 7, as expected. 
The only significant difference was that, for c* < 0, the electric field can change direction 
at t r = ~7/2a as mentioned above. If the reversal happens before the neutral line appears 
(t r < t c ) the final state is characterized by lower temperatures and streaming velocities 
than if the reversal happens after the appearance of the neutral line ( t r > t c ). 

Conclusions 

We see that the simple current sheet system we study here, which models the magne- 
totail fields by a Harris sheet with a small normal magnetic field component and a tearing 
perturbation added, and only considers the response of charged particles to these fields, 
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with no collective effects present other than the ones implicit in the use of the tearing 
mode, appears to describe correctly several features of the magnetotail plasma. Particle 
temperatures and bulk earthward drift velocities, both in the plasma sheet and in the 
PSBL, are in good agreement with the observed values, and the neutral line is seen to 
play an important role in the development of the tearing mode. The system also seems to 
always end up in the neighborhood of the same state shortly before the time the neutral 
line appears, and the further from that state it is when we start the simulation, the faster 
it will converge to it. This state appears to be characterized by /c~l. For runs that start 
with k 1, the value of k is drastically reduced in a very short time ( one particular ran 
went from k w 40 to k « 0.6 in less than 0.05 Toi = 0.16 sec) and then remains more or less 
unchanged, with k~1. Since k = by/ajpo ~ FJ -1 / 4 , this is consistent with the evolution of 
the temperature sis was discussed in Section-Ill. If we start with a higher temperature, so 
that the initial k is less than one, there is less heating sind the final temperature is again 
within the observed values, with k~1. Finally we have shown that orbitsd stochasticity can 
lead to the brosidening of the particle distribution function and hence irreversible heating, 
sind that a time dependent tearing mode amplitude has no effect on the results, particularly 
on the importance of the neutral line in the evolution of the magnetotail. 

The present work is confined to a zero cross-tsiil magnetic field component B y . This is 
not always true in the magnetotail, and there are indications (Tajima, 1981) that a nonzero 
By may inhibit the tearing mode, by magnetizing the electrons. This is clearly an important 
point that needs to be further clarified. The present work also does not address the trigger 
mechanism, nor can it properly treat the effect of the generated currents back onto the 
fields, although it suggests that neither a psirtial depletion, nor a reinforcement of the mode 
by the particle currents will introduce a qualitative change to our results. This is still an 
important question however, since the observed ion current densities are of the order of 
several fiA/m 2 , much greater than the background current density, jb ~ 2.5 x 10 -3 fiA/m 2 


20 


(although confined to a much smaller current channel). The observed electron current 
densities however are much smaller, of the order of 10 -5 fiA/m 2 . This suggest that the 
resolution of the point may be related to the work of Lembege and Pellat (1982), who 
showed that the adiabaticity of the electrons will also stabilize the ion tearing mode. If 
that work can be extended to show that the nonadiabaticity of the electrons can actually 
be viewed as a ‘valve’ that regulates the flow of free energy in the system, the ion current 
densities we observe will no longer present a problem. 

The fact that a test particle code with a simple model for the electromagnetic fields of 
the geomagnetic tail give such realistic results, and the fact that these results axe insensitive 
to initial conditions, or even to a small amount of velocity-space diffusion, lends support 
to the proposition put forth recently (Coroniti, 1985a, b; Buchner and Zelenyi, 1986, 1987; 
Martin, 1986) that stochastic particle dynamics in the magnetotail is an important ingre- 
dient of the tail behaviour. While self-consistent field effects are important in determining 
the global form of the electromagnetic perturbation of the geomagnetic tail as given by 
resistive MHD simulations, the rate of dissipation of the released magnetic energy and the 
rate of particle heating depend on the single particle processes discused here. 
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Figure Captions 


Figure- 1 

The magnetic field with 200 sample particles plotted for four different times for the 
ion run. (a) t=0 T 0 ,, (b) t=183.3 T oi , (c) t=233.3T oi , (d) t=280.0T oi 

Figure- 2 

(j • E) and temperature plots for the ion run. (a, b) t=166.7 To,; (c, d) t=216.7 Tor, 
(e, f) t=273.3 Tor (g, h) t=166.7 To,-, {j • E) and temperature plotted against x. 

Figure-3 

The earthward drift velocity v x for the ion rim, plotted for four different times, (a) 
t=33.3 T oi , (b) t=166.7 T 0 „ (c) t=246.7 T 0 „ (d) t=280.0 T 0 , 

Figure- 4 

Contour plots of the particle distribution function for the five boxed areas of particle 
concetration in Fig. (3)-c. a)-d), correspond to the four peripheral areas of particle conce- 
tration, clockwise from the top left, e) corresponds to the central area, —0.5 < z/a < 0.5, 
-0.2 < x/X < 0.2. 

Figure-5 

{j • E) plotted against time for the electron run. (a) The early run, no velocity-space 
diffusion, (b) The late run, no velocity-space diffusion, (c) The early run, D v = 0.6 km 2 /s 3 . 
(d) The late rim, D v = 0 .6 km 2 /s 3 . 
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